Method and system for processing an image featuring multiple scales

ABSTRACT

A method of processing an image is disclosed. The method comprises obtaining an image decomposed into a set of scaled images, each being characterized by a different image-scale; and calculating, for each of at least some scaled images, a relative luminance between the scaled image and another scaled image of the set, using intensities in the scaled image and intensities in the another scaled image. The method further comprises processing each scaled image using an adaptation procedure featuring an image-specific effective saturation function of the relative luminance, thereby providing a processed scaled image; combining at least some of the processed scaled images to provide a combined image; and outputting the combined image to a computer readable medium.

RELATED APPLICATIONS

This application is a Continuation-in-Part of PCT Patent Application No.PCT/IL2011/000639 having International filing date of Aug. 4, 2011,which claims the benefit of priority under 35 USC §119(e) of U.S.Provisional Patent Application No. 61/370,812 filed on Aug. 5, 2010. Thecontents of the above applications are all incorporated by reference asif fully set forth herein in their entirety.

FIELD AND BACKGROUND OF THE INVENTION

The present invention, in some embodiments thereof, relates to imageprocessing. Specifically, the present embodiments can be used forproviding an automatic dynamic range modulation of a digital image. Invarious exemplary embodiments of the invention the method and/orapparatus is used for companding (compressing and expanding) a highdynamic range (HDR) image. The present embodiments further comprise animaging system.

High dynamic range imaging (HDRI) is a set of techniques that allow afar greater dynamic range of exposures (large difference between lightand dark areas) than normal digital imaging techniques. The intention ofHDRI is to accurately represent the wide range of intensity levels foundin real scenes ranging from direct sunlight to the deepest shadows.

HDRI was originally developed for use with purely computer-generatedimages. Later, methods were developed to produce a HDR image from a setof photos taken with a range of exposures. With the rising popularity ofdigital cameras and easy to use desktop software, many amateurphotographers have used HDRI methods to create photos of scenes with ahigh dynamic range.

HDR images require a higher number of bits per color channel thantraditional images, both because of the linear encoding and because theyneed to represent values from 10⁻⁴ to 10⁸ (the range of visibleluminance values) or more. 16-bit (“half precision”) or 32-bit floatingpoint numbers are often used to represent HDR pixels. However, when theappropriate transfer function is used, HDR pixels for some applicationscan be represented with as few as 10-12 bits for luminance and 8 bitsfor chrominance without introducing any visible quantization artifacts.

Digital images may contain a huge amount of data, especially for highquality display and printing. Commercially available digital imagingdevices are known to acquire image information across a wide dynamicrange of several orders of magnitude. Additionally, there are softwaresolutions which fuse multiple exposures of the same scene at lowerdynamic range into one image of higher dynamic range.

Typically, although at the time of image capture the acquired dynamicrange is rather large, a substantial portion of it is lost once theimage is digitized, printed or displayed. For example, most images aredigitized to 8-bits (256 levels) per color-band, i.e., a dynamic rangeof about two orders of magnitude. The problem is aggravated once theimage is transferred to a display or a print medium which is oftenlimited to about 50 levels per color-band.

International Publication No. WO2009/081394, the contents of which arehereby incorporated by reference discloses an image processing techniquein which a digital HDR image is processed using two adaptationprocedures employed on the achromatic channel of the digital image. Eachadaptation procedure incorporates a different effective saturationfunction of the intensity. The adaptation procedures mimic the operationof the physiological visual system, wherein the first procedure mimicsthe “on” retinal pathway and the second adaptation procedure mimics the“off” retinal pathways. The intensity level of each picture-element ofthe digital image is processed by both procedures. The result of eachprocessing is an intermediate intensity level. All the intermediateintensity levels of the picture-element are then combined to provide anew achromatic intensity.

SUMMARY OF THE INVENTION

According to an aspect of some embodiments of the present inventionthere is provided a method of processing an image. The method comprises:obtaining an image decomposed into a set of scaled images, each beingcharacterized by a different image-scale; processing each scaled imageof the set using an adaptation procedure featuring an image-specificeffective saturation function of intensities in the scaled image andintensities in another scaled image of the set, thereby providing aprocessed scaled image; combining at least some of the processed scaledimages to provide a combined image; and outputting the combined image toa computer readable medium.

According to some embodiments of the present invention the methodcomprises, for each of at least some scaled image of the set,calculating a relative luminance between the scaled image and anotherscaled image of the set, using intensities in the scaled image andintensities in the another scaled image. According to some embodimentsof the present invention the image-specific effective saturationfunction is a function of the relative luminance.

According to some embodiments of the present invention the methodcomprises receiving the image and decomposing the image into the set ofscaled images.

According to some embodiments of the invention the decomposing comprisesselecting a size of the set based on a size of the image.

According to some embodiments of the invention the decomposing comprisesdetermining an amount of information in each scaled image being formed,and ceasing the decomposing when the amount of information is below apredetermined threshold.

According to an aspect of some embodiments of the present inventionthere is provided a method of capturing and displaying an image. Themethod comprises capturing an image of a scene and processing the imageusing the method as delineated above.

According to some embodiments of the present invention the methodcomprises recording radiation selected from the group consisting ofvisible light, infrared light, ultraviolet light, X-ray radiation,radiofrequency radiation, microwave radiation and ultrasound radiation,thereby capturing the image.

According to an aspect of some embodiments of the present inventionthere is provided a computer software product, comprising acomputer-readable medium in which program instructions are stored, whichinstructions, when read by a computer, cause the computer to execute themethod as delineated above.

According to an aspect of some embodiments of the present inventionthere is provided a system for processing an image, the system comprisesa data processor configured for: decomposing the image into a set ofscaled images, each being characterized by a different image-scale;processing each scaled image of the set using an adaptation procedurefeaturing an image-specific effective saturation function of intensitiesin the scaled image and intensities in another scaled image of the set,thereby providing a processed scaled image; and combining at least someof the processed scaled images to provide a combined image.

According to some embodiments of the present invention the dataprocessor calculates, for each scaled image of the set, a relativeluminance between the scaled image and another scaled image of the setusing intensities in the scaled image and intensities in the anotherscaled image. According to some embodiments of the present invention theimage-specific effective saturation function is a function of therelative luminance.

According to an aspect of some embodiments of the present inventionthere is provided an imaging system. The imaging system comprises animage capturing system and the processing system as delineated above.

According to some embodiments of the present invention the imagecapturing system is selected from the group consisting of a digitalcamera, a video camera, a CMOS digital camera, an infrared camera, anX-ray camera, a scanner, a microwave imaging, a computerized tomographyscanner, a magnetic resonance imaging scanner, a mammography scanner, anultrasonic scanner, an impedance imaging system, an endoscopic imagingdevice, a radio telescope, a digital telescope, a digital microscope anda system for translating an analog image to a digital image.

According to some embodiments of the present invention a characteristicdynamic range of the combined image is lower than a characteristicdynamic range of the original image.

According to some embodiments of the present invention the scaled imagesare combined by multiplication.

According to some embodiments of the present invention the set is anordered set and wherein the relative luminance is expressed as functionof a ratio between the intensities in the scaled image and theintensities in the other scaled image.

According to some embodiments of the present invention theimage-specific effective saturation function comprises an image-specificexponent, which is a function of a local contrast within thescale-image.

According to some embodiments of the present invention the processingcomprises modulating each relative luminance to provide a plurality ofmodulated relative luminance levels, wherein the combining comprisescombining the modulated relative luminance levels.

According to some embodiments of the present invention the modulatingcomprises selecting a relative luminance level such that two effectivesaturation functions corresponding to different image-specific exponentbut the same scale are substantially matched, more preferablysubstantially equal.

According to some embodiments of the present invention the localcontrast is calculated using a contrast-based adaptation procedureemployed for each picture-element of the scaled image.

According to some embodiments of the present invention thecontrast-based adaptation procedure calculates the local contrast basedon a difference between a second order opponent receptive field functioncalculated for the picture-element and a second order opponent receptivefield function calculated for nearby picture-elements.

According to some embodiments of the present invention theimage-specific exponent is a decreasing function of the local contrast.

According to some embodiments of the present invention theimage-specific exponent is a linear decreasing function of the localcontrast.

According to some embodiments of the present invention theimage-specific effective saturation function comprises a modulationfunction which is calculated based on a local contrast.

According to some embodiments of the present invention the modulationfunction has higher values when the local contrast is low, and lowervalues when the local contrast is high.

According to some embodiments of the present invention the methodcomprises employing a global gain operation for all scaled images of theset.

According to some embodiments of the present invention the global gainoperation features a global gain exponent, and the method comprisescalculating the global gain exponent using an optimization procedure.

According to some embodiments of the present invention the optimizationprocedure comprises selecting a set of candidate gain exponents,assigning a score to each candidate gain exponent, and selecting thegain exponent responsively to the score.

According to some embodiments of the present invention the scorecomprises a characteristic contrast.

According to some embodiments of the present invention the set is anordered set and wherein the scaled image and the other scaled image areadjacent images in the set.

According to some embodiments of the present invention the image is anHDR image.

According to some embodiments of the present invention the image is ofat least one type selected from the group consisting of a visible lightimage, a stills image, a video image, an X-ray image, an infrared image,a thermal image, a ultraviolet image, a computerized tomography (CT)image, a mammography image, a Roentgen image, a positron emissiontomography (PET) image, a magnetic resonance image, an ultrasoundimages, an impedance image, an elastography image, and a single photonemission computed tomography (SPECT) image.

Unless otherwise defined, all technical and/or scientific terms usedherein have the same meaning as commonly understood by one of ordinaryskill in the art to which the invention pertains. Although methods andmaterials similar or equivalent to those described herein can be used inthe practice or testing of embodiments of the invention, exemplarymethods and/or materials are described below. In case of conflict, thepatent specification, including definitions, will control. In addition,the materials, methods, and examples are illustrative only and are notintended to be necessarily limiting.

Implementation of the method and/or system of embodiments of theinvention can involve performing or completing selected tasks manually,automatically, or a combination thereof. Moreover, according to actualinstrumentation and equipment of embodiments of the method and/or systemof the invention, several selected tasks could be implemented byhardware, by software or by firmware or by a combination thereof usingan operating system.

For example, hardware for performing selected tasks according toembodiments of the invention could be implemented as a chip or acircuit. As software, selected tasks according to embodiments of theinvention could be implemented as a plurality of software instructionsbeing executed by a computer using any suitable operating system. In anexemplary embodiment of the invention, one or more tasks according toexemplary embodiments of method and/or system as described herein areperformed by a data processor, such as a computing platform forexecuting a plurality of instructions. Optionally, the data processorincludes a volatile memory for storing instructions and/or data and/or anon-volatile storage, for example, a magnetic hard-disk and/or removablemedia, for storing instructions and/or data. Optionally, a networkconnection is provided as well. A display and/or a user input devicesuch as a keyboard or mouse are optionally provided as well.

BRIEF DESCRIPTION OF THE DRAWINGS

The patent or application file contains at least one drawing executed incolor. Copies of this patent or patent application publication withcolor drawing(s) will be provided by the Office upon request and paymentof the necessary fee.

Some embodiments of the invention are herein described, by way ofexample only, with reference to the accompanying drawings and images.With specific reference now to the drawings in detail, it is stressedthat the particulars shown are by way of example and for purposes ofillustrative discussion of embodiments of the invention. In this regard,the description taken with the drawings makes apparent to those skilledin the art how embodiments of the invention may be practiced.

In the drawings:

FIG. 1 is a flowchart diagram illustrating a method suitable forprocessing an image, according to some embodiments of the presentinvention;

FIGS. 2A and 2B are plots of processed intensities R as a function of arelative luminance q, according to some embodiments of the presentinvention;

FIGS. 3A and 3B are schematic illustrations of rectangular grids ofpicture-elements which exemplify a concept of picture-element regions,according to various exemplary embodiments of the invention;

FIG. 4 is a schematic illustration of a system for processing an image,according to some embodiments of the present invention;

FIG. 5 is a schematic illustration of an imaging system, according tosome embodiments of the present invention;

FIGS. 6A and 6B show a thermal image before (6A) and after (6B)processing according to some embodiments of the present invention;

FIGS. 7A and 7B show another thermal image before (7A) and after (7B)processing according to some embodiments of the present invention;

FIGS. 8A and 8B show another thermal image before (8A) and after (8B)processing according to some embodiments of the present invention;

FIGS. 9A and 9B show another thermal image before (9A) and after (9B)processing according to some embodiments of the present invention;

FIGS. 10A and 10B show another thermal image before (10A) and after(10B) processing according to some embodiments of the present invention;

FIGS. 11A and 11B show another thermal image before (11A) and after(11B) processing according to some embodiments of the present invention;

FIGS. 12A-D show high dynamic range images, processed according to someembodiments of the present invention;

FIGS. 13A-F show a Gaussian pyramid obtained according to someembodiments of the present invention for a set of six scaled images;

FIGS. 14A-E show a luminance pyramid obtained according to someembodiments of the present invention from the Gaussian pyramid of FIGS.13A-F;

FIGS. 15A-E show a local contrast pyramid obtained according to someembodiments of the present invention from the luminance pyramid of FIGS.14A-E;

FIGS. 16A-E show a representation pyramid corresponding to a set ofexponents obtained according to some embodiments of the presentinvention from the local contrast pyramid of FIGS. 15A-E;

FIGS. 17A-E show a saturation pyramid obtained according to someembodiments of the present invention from the exponents corresponding tothe pyramid of FIGS. 16A-E; and

FIGS. 18A and 18B show an input image and an image processed accordingto some embodiments of the present invention using the pyramids of FIGS.13A-17E.

DESCRIPTION OF SPECIFIC EMBODIMENTS OF THE INVENTION

The present invention, in some embodiments thereof, relates to imageprocessing. Specifically, the present embodiments can be used forproviding an automatic dynamic range modulation of a digital image. Invarious exemplary embodiments of the invention the method and/orapparatus is used for companding (compressing and expanding) a highdynamic range (HDR) image. The present embodiments further comprise animaging system.

Before explaining at least one embodiment of the invention in detail, itis to be understood that the invention is not necessarily limited in itsapplication to the details of construction and the arrangement of thecomponents and/or methods set forth in the following description and/orillustrated in the drawings and/or the Examples. The invention iscapable of other embodiments or of being practiced or carried out invarious ways.

The present embodiments are concerned with method and system forprocessing an image to facilitate its display. At least part of theprocessing can be implemented by a data processing system, e.g., adedicated circuitry or a general purpose computer, configured forreceiving the image and executing the operations described below.

The method of the present embodiments can be embodied in many forms. Forexample, it can be embodied in on a tangible medium such as a computerfor performing the method operations. It can be embodied on a computerreadable medium, comprising computer readable instructions for carryingout the method operations. In can also be embodied in electronic devicehaving digital computer capabilities arranged to run the computerprogram on the tangible medium or execute the instruction on a computerreadable medium.

Computer programs implementing the method of the present embodiments cancommonly be distributed to users on a distribution medium such as, butnot limited to, a floppy disk, a CD-ROM, a flash memory device and aportable hard drive. From the distribution medium, the computer programscan be copied to a hard disk or a similar intermediate storage medium.The computer programs can be run by loading the computer instructionseither from their distribution medium or their intermediate storagemedium into the execution memory of the computer, configuring thecomputer to act in accordance with the method of this invention. Allthese operations are well-known to those skilled in the art of computersystems.

The image to be analyzed using the teachings of the present embodimentsis generally in the form of imagery data arranged gridwise in aplurality of picture-elements (e.g., pixels, group of pixels, etc.).

The term “pixel” is sometimes abbreviated herein to indicate apicture-element. However, this is not intended to limit the meaning ofthe term “picture-element” which refers to a unit of the composition ofan image.

References to an “image” herein are, inter alia, references to values atpicture-elements treated collectively as an array. Thus, the term“image” as used herein also encompasses a mathematical object which doesnot necessarily correspond to a physical object. The original andprocessed images certainly do correspond to physical objects which arethe scene from which the imaging data are acquired.

Each pixel in the image can be associated with a single digitalintensity value, in which case the image is a grayscale image.Alternatively, each pixel is associated with three or more digitalintensity values sampling the amount of light at three or more differentcolor channels (e.g., red, green and blue) in which case the image is acolor image. Also contemplated are images in which each pixel isassociated with a mantissa for each color channels and a common exponent(e.g., the so-called RGBE format). Such images are known as “highdynamic range” images.

The input image can be provided by any imaging modality, including,without limitation, a digital camera, a video camera, a CMOS digitalcamera, an infrared camera, a thermography device, an X-ray camera, ascanner, a microwave imaging, a computerized tomography scanner, asingle photon emission computed tomography device, a magnetic resonanceimaging scanner, a mammography scanner, an ultrasonic scanner, animpedance imaging system, an endoscopic imaging device, an elastographydevice, a radio telescope, a digital telescope, a digital microscope anda system for translating an analog image to a digital image.

Commercially available digital imaging devices based upon CCD detectorarrays are known to acquire image information across a wide dynamicrange of the order of 2 to 3 orders of magnitude. It is expected thatwith the rapid technologically development in the field of digitalimaging, this range will most likely be broadened in the near future.Typically however, although at the time of image capture the acquireddynamic range is rather large, a substantial portion of it is lost oncethe image is digitized, printed or displayed. For example, most imagesare digitized to 8-bits (256 levels) per color-band, i.e., a dynamicrange of about two orders of magnitude. The problem is aggravated oncethe image is transferred to a display or a print medium which is oftenlimited to about 50 levels per color-band.

A novel imaging technology, recently developed, employs CMOS with activepixel sensors [O. Yadid-Pecht and E. Fossum, “Image Sensor WithUltra-High-Linear-Dynamic Range Utilizing Dual Output CMOS Active PixelSensors”, IEEE Trans. Elec. Dev., Special issue on solid state imagesensors, Vol. 44, No. 10, 1721-1724], which are capable of locallyadjusting the dynamical range, hence to provide a high quality imagewith high dynamic range.

In addition, over the past years software solutions were developed forfuse multiple exposures of the same scene at low dynamic range (e.g.,256 levels per color-band) into one high dynamic range image (of about 4orders of magnitudes). High dynamic range images are typically providedin an RGBE format. In this format, 4 bytes are used (as opposed to 3bytes in conventional images) to create a representation similar tofloating point, where the first three bytes represent the three RGBcolor channels and the forth byte represents a common exponent to thethree colors channels. The dynamic range of such images is about 4orders of magnitude.

The motivation for developing imaging devices capable of capturing highdynamic range images is explained by the enormous gap between theperformances of the presently available devices and the ability of thehuman visual system to acquire detailed information from an ultra-highdynamic range scene. Specifically, the human visual system, which iscapable of acquiring a dynamic range of 14 orders of magnitude, caneasily recognize objects in natural light having a dynamic range of 12orders of magnitude.

Still, there is a growing gap between the state-of-the-art imagingdevices and display devices. High quality images, obtained either withphotographical film or by digital cameras, suffer, once displayed on ascreen or printed as a hard copy from loss in clarity of details andcolors at extreme light intensities, within shadows, dark regions,extremely bright regions and/or surfaces close to a lightening source.For example, as a single sharp edge in natural scene (e.g., a shadedobject in illuminated scene) can reach a dynamic range of 2 orders ofmagnitudes, presently available display devices may not be able torecognize such an edge. Another severe problem is that in a specificexposure a dark region of the image may be seen while a bright region isover exposed, or vise versa.

The technique developed by the present inventor is suitable for HDRimages as well as other images.

Referring now to the drawings, FIG. 1 is a flowchart diagramillustrating a method suitable for processing an image, according tosome embodiments of the present invention. The method of the presentembodiments can be used for processing any image including, withoutlimitation, a visible light image, a stills image, a video image, anX-ray image, a thermal image, a ultraviolet image, a computerizedtomography (CT) image, a mammography image, a Roentgen image, a positronemission tomography (PET) image, a magnetic resonance image, anultrasound images, an impedance image, and a single photon emissioncomputed tomography (SPECT) image.

The method begins at 10 and optionally continues to 11 at an image of ascene is captured. The image can be captured using any imaging techniqueknown in the art, including, without limitation, visible light imaging,infrared light imaging, ultraviolet light imaging, X-ray imaging,radiofrequency imaging, microwave imaging and ultrasound imaging. Theimaged scene can be of any type, including, without limitation, anoutdoor scene, an indoor scene, a nearby scene, a remote scene, anastronomical scene, an underwater scene, an intracorporeal scene (namelya scene that includes internal organs of a subject), an extracorporealscene (namely a scene that includes external organs of a subject), andany combination thereof.

Alternatively, 11 can be skipped in which case an image is received as astream of imaging data, as further detailed hereinabove.

In some embodiments of the present invention, the imaged is subjected toa preprocessing operation, such as, but not limited to, thepreprocessing operation described in International Patent ApplicationNo. PCT/IL2008/001419, the contents of which are hereby incorporated byreference. This embodiment is particularly useful when the image isformed by a computerized tomography technique, e.g., CT of SPECT.

At 12 the image is optionally and preferably decomposed into a set ofscaled images, each being characterized by a different image-scale.Alternatively, 12 can be skipped, in which case the set of scaled imagesis received by the method from an external source.

In various exemplary embodiments of the invention the set is an orderedset, wherein the kth element of the set is a blurred version of the k−1element. In other words, the images in the set are ordered such that theresolution of the k−1 image is finer than the resolution of the kthimage. The decomposition 12 can be done using any procedure known in theart.

A known operator for decomposing an image is referred to in theliterature as “Reduce.” A representative example of a Reduce operatorsuitable for the preset embodiments is described in Burt and Adelson,1983, “The Laplacian Pyramid as a Compact Image Code,” IEEE TransactionsOn Communications, vol. Com-31, No. 4, the contents of which are herebyincorporated by reference. When a Reduce operator is employed, the k+1element of the set can be calculated based on the kth element, asfollows I^(k+1)=Reduce(I^(k)). In some embodiments of the presentinvention the set of scaled images form an “Image PyramidRepresentation” which is formed by successively applying a downsamplingfilter having a weighting function centered on the pixel itself. In someembodiments of the present invention the weighting function is aunimodal function, such as a function having a shape which isapproximately a Gaussian. When the weighting function approximates aGaussian, the Image Pyramid Representation is referred to as a GaussianPyramid. However, it is not intended to limit the scope of the presentembodiments to unimodal or Gaussian weighting function. Other types ofweighting function, such as, but not limited to, a triangular weightingfunction and a trimodal weighting function, are not excluded from thescope of the present invention.

In some embodiments of the present invention a scaled image is obtainedby downsampling an image of a finer resolution. Thus, denoting theintensities of the kth scaled image by I^(k) (x, y), where the set oftuples (x, y) represents the picture-elements in the image, theintensities of the k+1 scaled image can be written as I^(k+1)(x,y)=I^(k)(ρ_(DS)x, ρ_(DS)y), where ρ_(DS) is a predetermined downsamplingcoefficient, and where I¹(x, y) can be, for example, the original image,denoted I_(in)(x,y). In some embodiments, the decomposing is done byintegrating the image I_(in)(x,y) with a kernel function, using adifferent spatial support for each resolution. A representative exampleof a kernel function for the kth scaled image is a Gaussian kernel,

${\frac{1}{\sqrt{\pi}\rho_{k}}{\exp\left( {- \frac{x^{2} + y^{2}}{\left( \rho^{k} \right)^{2}}} \right)}},$where ρ is a predetermined parameter.

In some embodiments, the decomposing employs an edge-preservingsmoothing filter, such as, but not limited to, a bilateraledge-preserving smoothing filter. A bilateral is a non-linear filterintroduced by Tomasi and Manduchi (see “Bilateral filtering for gray andcolor images,” Proc. IEE Intl. Conf. on Computer Vision, Bombay, India,1998), which is used for selective de-noising an images without blurringits edge. The bilateral filter takes into consideration both geometricdistances in the spatial domain and similarities in the intensitydomain. A bilateral filter typically features a convolution mask havingweights which are modified as a function of intensity differencesbetween a picture-element under consideration and its neighbors.

Denoting the edge-preserving smoothing filter by EPF, the kth scaledimage of the present embodiments can have the form EPF(I^(k−1)). Arepresentative example of a edge-preserving smoothing filter suitablefor the present embodiments is found, for example, in U.S. Pat. Nos.5,771,318, 7,146,059 and 7,199,793 the contents of which are herebyincorporated by reference.

As a representative example, which is not intended to be considered aslimiting, the following expression can be used as a bilateral filter:

$\begin{matrix}{{E\; P\;{F\left( {I\left( \overset{->}{r} \right)} \right)}} = \frac{\int{\int_{image}{{G_{r}\left( {{\overset{->}{r} - {\overset{->}{r}}^{\prime}}} \right)}{G_{s}\left( {{{I\left( {\overset{->}{r}}^{\prime} \right)} - {I\left( \overset{->}{r} \right)}}} \right)}{I\left( {\overset{->}{r}}^{\prime} \right)}{\mathbb{d}^{2}{\overset{->}{r}}^{\prime}}}}}{\int{\int_{image}{{G_{r}\left( {{\overset{->}{r} - {\overset{->}{r}}^{\prime}}} \right)}{G_{s}\left( {{{I\left( {\overset{->}{r}}^{\prime} \right)} - {I\left( \overset{->}{r} \right)}}} \right)}{\mathbb{d}^{2}{\overset{->}{r}}^{\prime}}}}}} & \left( {{EQ}.\mspace{14mu} 1} \right)\end{matrix}$where the vectors {right arrow over (r)} and {right arrow over (r)}′represent the coordinates of the picture-elements in the image, forexample, when the image is defined over a Cartesian coordinate system,{right arrow over (r)}=(x, y) and {right arrow over (r)}′=(x′, y′);G_(r) and G_(s) are localized functions with finite supports; and theintegration measure d²{right arrow over (r)}′ includes a Jacobian whichcorresponds to the coordinate system over which the image is defined,for example, for a Cartesian coordinate system d²{right arrow over(r)}′=dx′dy′. While some of the embodiments below are described, forclarity of presentation, by means of a Cartesian coordinate system, itis to be understood that more detailed reference to a Cartesiancoordinate system is not to be interpreted as limiting the scope of theinvention in any way. Notice that EQ. 1 features localized functionsboth in the spatial domain and in the intensity domain. Specifically,the function G_(r) is centered at the coordinate {right arrow over (r)}and the function Gs is centered at the intensity I({right arrow over(r)}) that is associate with the coordinate {right arrow over (r)}.

The localized functions G_(r) and G_(s) can have any form provided ithas a finite support. Representative examples including, withoutlimitation, Gaussians, Lorenzians, modified Bessel functions. In someembodiments of the present invention both G_(r) and G_(s) are Gaussians,e.g.,G _(s)(|r|)=exp(−r ²σ_(s) ²);(G _(r)(r)=exp(−r ²/σ_(r) ²),  (EQ. 2)where σ_(s) and σ_(r) are radius parameters characterizing the localsupport of G_(s) and G_(r), respectively.

The decomposition 12 preferably features a set of image-specificfilters, wherein for each scaled image an image-specific filter withdifferent image-specific parameters is employed. For example, when EQs.1 and 2 are employed, each scaled image is associated with a differentset of radius parameters. The radius parameters used for obtaining thekth image I^(k) are denoted σ_(r) ^(k) and σ_(s) ^(k). Thus, the kthimage t is preferably calculated using the expressionI^(k)=EPF(I^(k−1)), wherein EPF features two localized functions G_(r)and G_(s) with two respective image-specific radius parameters σ_(r)^(k) and σ_(s) ^(k).

The method optionally and preferably continues to 13 at which each of atleast some of the scaled images is processed using an adaptationprocedure featuring an image-specific effective saturation function.

The effective saturation function is “image-specific” in the sense thatfor each scaled image the procedures defined a specific effectivesaturation function which is typically different from the effectivesaturation function defined for any other scaled image in the set. Invarious exemplary embodiments of the invention the effective saturationfunction is applied for each picture-element of the scaled image beingprocessed and can therefore be viewed as a processed scaled image. Inother words, the returned values of the effective saturation functioncan be used as image intensities. The effective saturation function forthe kth scaled image I^(k)(x,y) is denoted R^(k)(x,y), and isinterchangeable referred to herein as the kth processed scaled image.

R^(k) is optionally and preferably a function of intensities in kthimage as well of intensities in at least one scaled image of the setwhich is other than the kth scaled image. In various exemplaryembodiments of the invention R^(k) is a function of intensities in atleast the kth image and an image which is a blurred version (e.g., witha coarser resolution) of the kth image. For example, R^(k) can be afunction of the I^(k) and I^(k+1). A typical expression for R^(k) is(for clarity of presentation, the spatial dependence of R^(k), I^(k) andI^(k+1) on the location (x, y) of the picture-element has been omitted):

$\begin{matrix}{R^{k} = \frac{R_{\max}I^{k}}{{\alpha\; I^{k}} + {\beta\; I^{k + 1}}}} & \left( {{EQ}.\mspace{14mu} 3} \right)\end{matrix}$where, R_(max), α and β are coefficients, which can be constants or theycan vary across the image and/or between images. For example, R_(max), αand β can each be set to 1, but other values are not excluded from thescope of the present invention.

In some embodiments, R^(k) is a function of a relative luminance q^(k)between the kth scaled image and the other scaled image (e.g., the k+1scaled image). In various exemplary embodiments of the invention q^(k)is the relative luminance per picture-element, namely it has a specificvalue for each picture-element in the kth scaled image. In theseembodiments, q^(k)=q^(k)(x,y), where the set of tuples (x, y) representsthe picture-elements in the kth image. The relative luminance q^(k)(x,y)can optionally and preferably be expressed in terms of the ratio betweenintensities in the kth scaled image and intensities in the other scaledimage. For example, q^(k)(x,y) can be defined asq^(k)(x,y)=f(I^(k)(x,y)/I^(k+1)(x,y)), where f is some function,preferably a monotonically increasing function, e.g., a linear functioncharacterized by a positive slope. Thus, in some embodiments of thepresent invention q^(k)(x,y) is defined as a I^(k)(x,y)/I^(k+1)(x,y)+b,where a and b are parameters which are optionally constants. In someembodiments, a=1 and b=0, but other values are not excluded from thescope of the present invention.

The calculation of the relative luminance q^(k) optionally andpreferably comprises some interpolation of the picture-elements in thecoarser image as known in the art. For example q^(k) can be calculatedas f(I^(k)(x,y)/Expand(I^(k+1)(x,y))), where Expand is an interpolationoperator that preferably matches the number of elements in I^(k+1) withthe number of elements in I^(k) using an interpolation algorithm. Insome embodiments of the present invention Expand is the reverse operatorof the Reduce operator. A representative example of an Expand operatorsuitable for the present embodiments is described in Burt and Adelson,supra.

A typical form of R^(k) when expressed as a function of the relativeluminance q^(k) is (for clarity of presentation, the spatial dependenceof q^(k) and R^(k) has been omitted):

$\begin{matrix}{{R^{k}\left( q^{k} \right)} = {\frac{R_{\max}}{\alpha + \left( {M/q^{k}} \right)^{\gamma}} + B}} & \left( {{EQ}.\mspace{14mu} 4} \right)\end{matrix}$where R_(max) and α are parameters already introduced above, M and γ area modulation coefficient and an exponent, respectively, and B is anoffset parameter. Each of M, γ and B can be a constant or some functionof the intensity. In the simplest case, α, M and γ are set to 1, and Bis set to zero, so that R^(k) is reduced to the formR^(k)=Rmax/(a+(q^(k))⁻¹). However, this need not necessarily be thecase, since, for some applications, it may be desired to let γ and/or Mand/or B be different from the above values and/or vary. Higher valuesof γ cause enhancement of the rate of change in R as a function of q,particularly at the vicinity of q=1, where there are small differencesbetween the luminance of a specific location and its context. The effectof the exponent γ on R^(k) is exemplified in FIG. 2A, which are plots ofR^(k) as a function of q^(k), for Rmax=M=1, B=0 and seven fixed valuesof γ: γ=0.5, γ=1, γ=2, γ=5, γ=7, γ=10 and γ=15. It is to be understoodthat these values are for illustrative purpose only and are not to beconsidered as limiting.

In some embodiments, the value of y is specific to the scaled image. Inthese embodiments, the exponent used for the kth scaled image is denotedγ^(k). In some embodiments of the present invention γ^(k) is adecreasing function, e.g., a linear decreasing function, of k.Preferably, γ^(k) is positive for all values of k.

In some embodiments, for the coarse scales (high k) that reflect theillumination γ^(k) can be less than 1 so as to compress the high dynamicrange of the illumination; at finer scales (low k) γ^(k) can be set to avalue which is higher than 1. Alternatively, γ^(k) is above 1 for allresolution. In some embodiments, γ^(k) satisfies γ_(min)≦γ^(k)≦γ_(max),where γ_(min) and γ_(max) are predetermined parameters which are thesame for all scales. In some specific embodiments of the presentinvention γ_(min) is from about 1 to about 3, e.g., about 2 and γ_(max)is from about 5 to about 9, e.g., about 7.

As a representative and non limiting example for a linear decrease of γas a function of the resolution index k, γ can be decreased by Δγ foreach integer increment of k, where Δγ is from about 0.1 to about 0.4 orfrom about 0.2 to about 0.3, e.g., about 0.25.

The modulation coefficient M can be viewed as a parameter whichmodulates the relative luminance q^(k). Formally, the expression q^(k)/Mcan be defined as an effective relative luminance, wherein higher valuesof M correspond to lower effective relative luminance and lower valuesof M correspond to higher effective relative luminance. The effect ofthe modulation coefficient M on R^(k) is exemplified in FIG. 2B, whichare plots of R^(k) as a function of q^(k), for R_(max)=γ=1 and threefixed values of M: M=0.5, M=1 and M=2. It is to be understood that thesevalues are for illustrative purpose only and are not to be considered aslimiting. Generally, larger values of M suppress the value of R^(k). Mcan be a global constant or it can vary over the picture-elements of thescaled image being processed and/or across the scaled images in the set.When M varies over the picture-elements it is realized as a function ofthe coordinates, for example, M=M(x,y), and when M varies across thescaled images of the set, the method features a set {M^(k)} ofcoefficients each of which can be a constant coefficient or a function,e.g., M^(k)=M^(k)(x,y).

While the embodiments above are described with a particular emphasis toan effective relative luminance having the form expression q^(k)/M, itis to be understood that more detailed reference to such expression isnot to be interpreted as limiting the scope of the invention in any way.Generally, an effective relative luminance {circumflex over (q)}^(k) canbe obtained using any linear or non-linear modulation operation, inwhich case the effective saturation function can be written as:

$\begin{matrix}{{R^{k}\left( q^{k} \right)} = {\frac{R_{\max}}{\alpha + \left( {\hat{q}}^{k} \right)^{- \gamma}} + B}} & \left( {{{EQ}.\mspace{14mu} 4}A} \right)\end{matrix}$

In some embodiments of the present invention the value of the exponent γand/or coefficient M is calculated based on image intensities in thescaled image.

For example, the image-specific exponent γ^(k) can be a function of alocal contrast C^(k) within the scaled image I^(k). Preferably, γ^(k)decreases with C^(k). In some embodiments, γ^(k) decreases linearly withC^(k), and in some embodiments γ^(k) decreases non-linearly with C^(k).A preferred relation between γ^(k) and C^(k) is:γ^(k) =f(C _(max))−C ^(k),  (EQ. 5)where C_(max) is a constant parameter which, in some embodiments, is themaximal contrast over image I^(k), and f is some function.Representative examples of expressions suitable for the functionf(C_(max)) including, without limitation, f(C_(max))=C_(max), andf(C_(max))=pC_(max), where p is a constant parameter which is preferablylarger than 1.

An alternative expression for a linear relation between γ^(k) and C^(k)is:γ^(k)=δ(1−C ^(k))  (EQ. 5A)where δ is a contrast enhancement parameter. Typically, larger values ofδ correspond to higher characteristic contrast of the final image.

Also contemplated is a non-linear relation between γ^(k) and C^(k), forexample,γ^(k) =N/(C ^(k))^(n),  (EQ. 6)where N and n are positive constants, e.g., N=n=1.

The modulation operation executed for providing the effective relativeluminance {circumflex over (q)}^(k) (for example, the value of M, inembodiments in which {circumflex over (q)}^(k)=q^(k)/M) is optionallyand preferably selected such that the returned values of the effectivesaturation functions for two different exponents be approximately equalfor a given scale. For example, consider the kth the effectivesaturation function R^(k). This function can be calculated more thanonce, using a different value for the exponent γ^(k) at eachcalculation. Without loss of generality, suppose that for a given scalek, two saturation functions R₀ ^(k) and R₁ ^(k) are calculated, asfollows:

${{R_{0}^{k}\left( q^{k} \right)} = {\frac{R_{\max}}{\alpha + \left( {M/q^{k}} \right)^{y_{0}^{k}}} + B}};$${R_{1}^{k}\left( q^{k} \right)} = {\frac{R_{\max}}{\alpha + \left( {M/q^{k}} \right)^{\gamma_{1}^{k}}} + {B.}}$

The exponents γ₀ ^(k) can have a fixed and predetermined (e.g., the samevalue for γ₀ ^(k) for all values of k), and γ₁ ^(k) can be selectedaccording to the local contrast as further detailed hereinabove. Invarious exemplary embodiments of the invention γ₀ ^(k) is the lowestallowed value of the exponent, e.g., γ₀ ^(k)=γ_(min). In someembodiments of the present invention the effective relative luminance{circumflex over (q)}^(k) is selected such that R₀ ^(k)({circumflex over(q)}^(k))=R₁ ^(k)(q^(k)).

The local contrast C^(k) can be calculated from the intensity values ofthe picture-element in the respective scaled image using any knownprocedure for detecting or calculating local contrast. Representativetechniques suitable for the present embodiments are found, for example,in U.S. Pat. Nos. 7,791,652, 7,929,739, 6,078,686 and 5,838,835 thecontents of which are hereby incorporated by reference.

In some embodiments of the present invention the local contrast iscalculated based on one or more intensity differences between scaledimages in the set. Typically, the local contrast is calculated based onintensity differences between scaled images whose resolution is notlower than the resolution of the currently processed scaled image. Arepresentative example for the local contrast C^(k) in these embodimentsis:

$\begin{matrix}{{C^{k} = {\sum\limits_{i = 0}^{k}{{I^{i} - I^{i + 1}}}^{ɛ}}},} & \left( {{{EQ}.\mspace{14mu} 6}A} \right)\end{matrix}$where ε is a local contrast exponent. Typical values for ε are fromabout 0.1 to about 1, e.g., about 0.3.

In some embodiments of the present invention the local contrast iscalculated using a contrast-based adaptation procedure which can beconstructed so as to mimic a mechanism of the human vision system knownas a second order achromatic induction. The contrast-based adaptationprocedure of the present embodiments is preferably as follows.

Firstly, the procedure mimics the transformation of a visual stimulusinto a response of post retinal second order opponent receptive fields(SORF's). These SORF's may refer to cortical levels even though thereceptive fields are not necessarily oriented. These SORF's have variousspatial resolutions, in compliance with the diversity found in the humanvisual system. The number of different spatial resolutions employed bythe procedure and the number of scaled images in the set can be the sameor they can be different. In some embodiments of the present inventionthe number of different spatial resolutions employed by thecontrast-based procedure is larger than the number of scaled images inthe set.

Secondly, local and remote contrasts are calculated based on the multiscale SORF responses, and thirdly a contrast-contrast induction isemployed. The contrast-contrast induction serves as a contrast gaincontrol and is expressed by the adapted responses of the SORF cells.

In the human visual system, the SORF cells receive their input from theretinal ganglion cells through several processing layers. The retinalganglion cells perform a (first order) adaptation and the SORF cellsreceive their responses after the adaptation. In the followingdescription, the first order adaptation is not modeled for clarity ofpresentation, but the skilled artisan, provided with the informationdescribed herein would know how to employ first order adaptation, e.g.,using the formalism of center and surround adaptation terms describedabove and/or the first order adaptation described in Barkan et al.,(2008), “Computational adaptation model and its predictions for colorinduction of first and second orders,”, J. of vision 8(7) 27 1-26.

The SORF cells have an opponent type receptive field with acenter-surround spatial structure. Thus, in various exemplaryembodiments of the invention the method defines, for eachpicture-element 20 one or more regions in the vicinity ofpicture-element 20 (but not necessarily adjacent thereto). Typically, asurrounding region is defined. In some embodiments the method alsodefines a center region which comprises element 20 and picture elementsimmediately adjacent to element 20. Alternatively, the center region maybe a single element region, hence comprising only element 20. Thisalternative, of course, coincide with the embodiment in which no centerregion is selected.

The concept of the center and surrounding regions may be betterunderstood from the following example, with reference to FIGS. 3A-B.Thus, if the picture elements are arranged in a rectangular grid 30, thecenter region may be a single picture element (element 20), and thesurround region may be picture elements 32 surrounding picture elements20. Picture elements 34, surrounding picture elements 32 can be referredto as remote region.

In FIG. 3A, the surround region comprises eight picture-elementsimmediately surrounding (i.e., adjacent to) element 20, and the remoteregion comprises 40 picture-element forming the two layers surroundingthose eight picture-elements. However, this need not necessarily be thecase, since, for some applications, it may be desired to extend thesurround region farther from those eight elements which immediatelysurround element 20. FIG. 3B, for example, illustrates an embodiment inwhich the surround region comprises 48 picture-element which form thefirst three layers surrounding element 20, and the remote regioncomprises 176 picture-element which form the four layers surroundingthose 48 elements. Also contemplated are embodiments in which the centerregion comprises more picture-element, e.g., an arrangement of 2×2 or3×3 picture-elements. Other definitions for the center, surrounding andremote regions are not excluded from the present invention, both for arectangular grid or for any other arrangement according to which thepicture elements of the scaled image are received.

Once the region(s) are defined, the intensities of the picture elementsin each region are preferably used for calculating, for each region, anoverall regional intensity. The overall intensity can be calculated as aconvolution of the intensity of the picture elements in each region withthe respective regional spatial profile. For the center region, thisconvolution is preferably realized by the following equations:

$\begin{matrix}{L_{cen}^{k} = {\int{\int_{cen}{{I^{k}\left( {x,y} \right)}{f_{c}^{k}\left( {{x - x_{0}},{y - y_{0}}} \right)}{\mathbb{d}x}{\mathbb{d}y}}}}} & \left( {{EQ}.\mspace{14mu} 7} \right)\end{matrix}$where (x₀, y₀) is the location of the center of the center region, f_(c)^(k) is the spatial profile for the center region at the kth scaledimage, and the integration is over all the picture-elements in thecenter region. Without loss of generality, (x₀, y₀) can be set to (0,0). In the following description, the reference to x₀ and y₀ istherefore omitted. The center spatial profile is preferably a spatialdecaying function. Examples for the specific form of the spatial profileof the center region include, but are not limited to, a Gaussian, anexponent, a Lorenzian, a modified Bessel function and a power-decayingfunction.

In some embodiments of the present invention, f_(c) ^(k) is given by:

$\begin{matrix}{{{{f_{c}^{k}\left( {x,y} \right)} = \frac{\exp\left( {{- x^{2}} + {y^{2}/\rho_{cen}^{k\; 2}}} \right)}{\pi \cdot \rho_{cen}^{k}}};x},{y \in {center\_ area}}} & \left( {{EQ}.\mspace{14mu} 8} \right)\end{matrix}$where ρ_(cen) represents the radius of the center region.

For the surround region, the convolution is similarly defined:

$\begin{matrix}{L_{srnd}^{k} = {\int{\int_{srnd}{{I^{k}\left( {x,y} \right)}{f_{s}^{k}\left( {x,y} \right)}{\mathbb{d}x}{\mathbb{d}y}}}}} & \left( {{EQ}.\mspace{14mu} 9} \right)\end{matrix}$where f_(c) ^(k) is the spatial profile for the surround region at thekth scaled image, and the integration is over all the picture-elementsin the surround region. Examples for the specific form of the spatialprofile of the surround region include, but are not limited to, aGaussian, an exponent, a Lorenzian, a modified Bessel function and apower-decaying function.

In some embodiments of the present invention, f_(s) ^(k) is given by:

$\begin{matrix}{{{{f_{s}^{k}\left( {x,y} \right)} = \frac{\exp\left( {{- x^{2}} + {y^{2}/\rho_{sur}^{k\; 2}}} \right)}{\pi \cdot \rho_{sur}^{k}}};x},{y \in {surround\_ area}}} & \left( {{EQ}.\mspace{14mu} 10} \right)\end{matrix}$where ρ_(sur) (also referred to below as ρ_(sur)) represents the radiusof the surrounding region. The total weight of f_(c) and f_(s) istypically 1.

The response of the SORF L_(sorf) can be calculated based on thedifference between the intensity value of picture-element(s) in thecenter region and intensity values of nearby picture-elements in thesurround region. For example,L _(sorf) ^(k) =L _(cen) ^(K) −L _(srnd) ^(k)  (EQ. 11)

The calculations of the SORF response can be alternatively derived froman achromatic double-opponent receptive field structure, where thecenter and the surrounding regions contained opponent sub-unitsdescribed for example, in Barkan et al. 2008, supra.

L^(k) _(sorf) is optionally and preferably used for calculating acontrast term, for example, by averaging the value of L^(k) _(sorf) overseveral scaled images or, more preferably, over all the scaled images inthe set. A representative examples for the contrast C is given by:

$\begin{matrix}{{C^{k}\left( {x,y} \right)} = {\sum\limits_{k^{\prime}}\frac{\int{\int{{{L_{sorf}^{k^{\prime}}\left( {x^{\prime},y^{\prime}} \right)}}^{\delta}{W^{k^{\prime}}\left( {{x^{\prime} - x},{y^{\prime} - y}} \right)}{\mathbb{d}x}{\mathbb{d}y}}}}{\int{\int{{W^{k^{\prime}}\left( {x,y} \right)}{\mathbb{d}x}{\mathbb{d}y}}}}}} & \left( {{EQ}.\mspace{14mu} 12} \right)\end{matrix}$where the integrations are preferably over a region encompassing thescales that are defined in the contrast region, δ is a parameter, andW^(k)(x,y) are weight functions which are preferably localized at (x,y),with some predetermined support. The δ can be any integer or non-integerpositive number. Typically, the δ parameter satisfies δ≧1, e.g., δ=1 orδ=2 or δ=3. A representative example of W^(k)(x,y) is a two-dimensionalGaussian:

$\begin{matrix}{{{W^{k}\left( {x,y} \right)} = {\exp\left( {- \frac{x^{2} + y^{2}}{\rho_{local}^{2}}} \right)}},} & \left( {{EQ}.\mspace{14mu} 13} \right)\end{matrix}$where ρ_(local) is a radius parameter representing the size of thesupport.

Once calculated, C^(k) can be used for calculating the exponent γ^(k),by means of a linear or non-linear relation, as further detailedhereinabove.

The above contrast-based procedure can also be used for calculating themodulation coefficient M described above (see EQ. 4). In someembodiments, M is a decreasing function of L^(k) _(sorf). In other wordsM has higher values when the local contrast is low, and lower valueswhen the local contrast is high. For example, the coefficient M^(k) forthe kth scaled image can be a linear decreasing function of L_(sorf),e.g., M=f(I^(k) _(max))−L^(K) _(sorf), where I^(k) _(max) is the maximalintensity over the kth scaled image, and f is some function, e.g.,f(I^(k) _(max))=I^(k) _(max) or f(I^(k) _(max))=mI^(K) _(max) where m isa positive parameter. Alternatively, M can be a non-linear decreasingfunction of L_(sorf). A representative example is the expressionM=1/L^(k) _(sorf) or the like. Other expressions are not excluded fromthe scope of the present invention.

When EQ. 3 is employed, the contrast-based procedure can be used forcalculating the β coefficient. In some embodiments, the β coefficient isa decreasing function of L^(k) _(sorf). In other words β has highervalues when the local contrast is low, and lower values when the localcontrast is high. For example, the coefficient β^(k) for the kth scaledimage can be a linear decreasing function of L_(sorf), e.g.,β^(k)=f(I^(k) _(max))−L^(k) _(sorf), where I^(k) _(max) is the maximalintensity over the kth scaled image, and f is some function, e.g.,f(I^(k) _(max))=I^(k) _(max) or f(I^(k) _(max))=mI^(k) _(max) where m isa positive parameter. Alternatively, M can be a non-linear decreasingfunction of L_(sorf). A representative example is the expressionβ^(k)=1/L^(k) _(sorf) or the like. Other expressions are not excludedfrom the scope of the present invention.

Referring again to FIG. 1 the method can continue to 14 at which atleast some of the processed scaled images are combined to provide acombined image. The images are preferably combined by multiplication.For example a combined image I_(combined) can be obtained using thefollowing equation:

$\begin{matrix}{{I_{combined} = {\prod\limits_{k}I^{k}}},} & \left( {{EQ}.\mspace{14mu} 14} \right)\end{matrix}$where the multiplication is over some or all the scaled images in theset.

The combined image can also be obtained from the relative luminance ofthe effective saturation functions, rather than from the effectivesaturation functions themselves. These embodiments are particularlyuseful when the relative luminance levels are modulated. For example,the combined image can be obtained by multiplying the effective relativeluminance values of at least a few, more preferably, all the effectivesaturation functions:

$\begin{matrix}{{I_{combined} = {\prod\limits_{k}{\hat{q}}^{k}}},} & \left( {{{EQ}.\mspace{14mu} 14}A} \right)\end{matrix}$where the multiplication is over some or all the scaled images in theset.

The calculation of the combined image optionally and preferablycomprises some interpolation of the picture-elements in the coarserimage as known in the art. For example, the following iterative processcan be employed:I ^(n) ={circumflex over (q)} ^(n) {circumflex over (q)} ^(n+1){circumflex over (q)} ^(n+1)=Expand(I ^(n+1)),  (EQ. 14B)where Expand is an interpolation operator as further detailedhereinabove. Once the iterative process (14B) is completed, the combinedimage can be defined as image obtained at the last iterative step, e.g.,I⁰.

The method optionally continues to 15 at which the combined image isnormalized. A normalization procedure suitable for the presentembodiments is the so called log-mean normalization, wherein thelogarithm of the intensity of each picture-element is first normalizedby the average of the logarithms of intensities and then exponentiated.Formally, this procedure can be written as:

$\begin{matrix}{{I_{combined}->{\exp\left\lbrack {{\log\left( I_{combined} \right)}\frac{TL}{\left\langle {\log\left( I_{combined} \right)} \right\rangle}} \right\rbrack}},} & \left( {{EQ}.\mspace{14mu} 15} \right)\end{matrix}$where <log(I_(combined))> is the average of the logarithms ofintensities calculated over the combined image, and TL is a constantparameter which represents the log-mean of the normalized image. Themethod optionally and preferably proceeds to 16 at which the combinedand/or normalized image is transmitted to a computer readable medium,from which it can be displayed or printed as desired.

In various exemplary embodiments of the invention the characteristicdynamic range of the combined and/or normalized image is lower than thecharacteristic dynamic range of the original image. For example, whenthe characteristic dynamic range of the original image spans over 4 ormore orders of magnitudes, the characteristic dynamic range of thecombined image is three or two orders of magnitudes. In variousexemplary embodiments of the invention the characteristic dynamic rangeof the combined image is sufficient to allow all the intensities of theimage to be displayed on a display device. For example, in someembodiments of the present invention the combined image comprises nomore than 256 different intensities.

The method ends at 17.

FIG. 4 is a schematic illustration of a system 40 for processing animage, according to some embodiments of the present invention. System 40comprises a data processor 42 having a computation module whichcomprises at least an image decomposer 44 configured for decomposing theimage into a set of scaled images, a scaled image processing module 46configured for processing each scaled image of set to provide aprocessed scaled image, and an image combiner 48 configured forcombining at least some of the processed scaled images to provide acombined image. In various exemplary embodiments of the invention thecomputation module of data processor 42 is configured for executing atleast some of the operations described above with respect to method10-17.

FIG. 5 is a schematic illustration of an imaging system 50, according tosome embodiments of the present invention. Imaging system 50 comprisesan image capturing system 52 and system 40. Image capturing system 52can be of any type, including, without limitation, a digital camera, avideo camera, a CMOS digital camera, an infrared camera, an X-raycamera, a scanner, a microwave imaging, a computerized tomographyscanner, a magnetic resonance imaging scanner, a mammography scanner, anultrasonic scanner, an impedance imaging system, an endoscopic imagingdevice, a radio telescope, a digital telescope, a digital microscope anda system for translating an analog image to a digital image.

As used herein the term “about” refers to ±10%.

The word “exemplary” is used herein to mean “serving as an example,instance or illustration.” Any embodiment described as “exemplary” isnot necessarily to be construed as preferred or advantageous over otherembodiments and/or to exclude the incorporation of features from otherembodiments.

The word “optionally” is used herein to mean “is provided in someembodiments and not provided in other embodiments.” Any particularembodiment of the invention may include a plurality of “optional”features unless such features conflict.

The terms “comprises”, “comprising”, “includes”, “including”, “having”and their conjugates mean “including but not limited to”.

The term “consisting of” means “including and limited to”.

The term “consisting essentially of” means that the composition, methodor structure may include additional ingredients, steps and/or parts, butonly if the additional ingredients, steps and/or parts do not materiallyalter the basic and novel characteristics of the claimed composition,method or structure.

As used herein, the singular form “a”, “an” and “the” include pluralreferences unless the context clearly dictates otherwise. For example,the term “a compound” or “at least one compound” may include a pluralityof compounds, including mixtures thereof.

Throughout this application, various embodiments of this invention maybe presented in a range format. It should be understood that thedescription in range format is merely for convenience and brevity andshould not be construed as an inflexible limitation on the scope of theinvention. Accordingly, the description of a range should be consideredto have specifically disclosed all the possible subranges as well asindividual numerical values within that range. For example, descriptionof a range such as from 1 to 6 should be considered to have specificallydisclosed subranges such as from 1 to 3, from 1 to 4, from 1 to 5, from2 to 4, from 2 to 6, from 3 to 6 etc., as well as individual numberswithin that range, for example, 1, 2, 3, 4, 5, and 6. This appliesregardless of the breadth of the range.

Whenever a numerical range is indicated herein, it is meant to includeany cited numeral (fractional or integral) within the indicated range.The phrases “ranging/ranges between” a first indicate number and asecond indicate number and “ranging/ranges from” a first indicate number“to” a second indicate number are used herein interchangeably and aremeant to include the first and second indicated numbers and all thefractional and integral numerals therebetween.

It is appreciated that certain features of the invention, which are, forclarity, described in the context of separate embodiments, may also beprovided in combination in a single embodiment. Conversely, variousfeatures of the invention, which are, for brevity, described in thecontext of a single embodiment, may also be provided separately or inany suitable subcombination or as suitable in any other describedembodiment of the invention. Certain features described in the contextof various embodiments are not to be considered essential features ofthose embodiments, unless the embodiment is inoperative without thoseelements.

Various embodiments and aspects of the present invention as delineatedhereinabove and as claimed in the claims section below find experimentalsupport in the following examples.

EXAMPLES

Reference is now made to the following examples, which together with theabove descriptions illustrate some embodiments of the invention in a nonlimiting fashion.

Example 1 Thermal Images

The technique of the present embodiments was employed to process thermalimages acquired using an infrared camera. FIGS. 6A, 7A, 8A, 9A, 10A and11A show the original images, before processing. Each image wasprocessed according to some embodiments of the present invention asdescribed in the flowchart of FIG. 1. The image was decomposed into 3scaled images each of which was processed using EQs. 1, 2, 4, 5 and7-14. The modulation coefficient was set to 1.

For the contrast-based adaptation procedure, up to six possibleresolutions were employed. The application of the specific resolutionwas done such that the coarse resolution was applied to the overallscales for the whole image. In other word, fine resolution in the SORFwere not included in the coarse general scales. Table 1 specifies theradii of the center and surround regions which were used for eachresolution k.

TABLE 1 k radius of the center region radius of the surround region 1 13 2 3 9 3 5 15 4 7 12 5 9 27 6 11 33

FIGS. 6B, 7B, 8B, 9B, 10B and 11B show the images after processing adcombining. As shown, the technique of the present embodiments allowsdistinguishing between image features which were undistinguishablebefore the processing.

Example 2 HDR Images

Embodiments of the present invention were applied to High Dynamic Rangeimages in an RGBE format. The original HDR images were obtained throughthe courtesy of Michael Werman, Erik Reinhard, Greg Ward, SpheronVR AG,Munsell Color Science Laboratory, and Paul Debevec's website.

Achromatic intensities were extracted from the polychromatic data of theimages. This was performed by transforming each pixel in the RGBE imageto CIE XYZ using the D65 sRGB transform matrix [IEC 61966-2-1:1999]. Theachromatic intensity of each pixel was defined as the Y value of thepixel. Further transformation from the CIE XYZ space to CIE xyz wasperformed for each pixel, and the x,z, values were applied to the newachromatic intensities yielded according to various exemplaryembodiments of the present invention.

Each image was processed according to some embodiments of the presentinvention as described in the flowchart of FIG. 1. The image wasdecomposed into 4 scaled images each of which was processed using EQs.1, 2, 4, 5 and 7-14. The modulation coefficient was set to 1. Theparameters used in the contrast-based adaptation procedure were asdetailed in Example 1, above.

The dynamic ranges of original images are larger than the maximaldisplayable dynamic range of a conventional display device or aconventional printer and are therefore not shown. The combined andnormalized processed scaled images are shown in FIGS. 12A-D.

The present example demonstrated the ability of embodiments of thepresent invention to perform automatic high dynamic range compression.The difference between the dynamic ranges of the original and processedimage was up to about 10¹⁰ levels of intensity. The results demonstratea significant compression while preserving, and even slightly enhancingthe details in the images, both in the bright and dark zones. Thetechnique has been applied for a large number of images. Although mostof the experiments were performed using the same set of parameters, thedynamic compression was successful in all processed images. Yet, thetechnique of the present embodiments can be applied to different extentsby assigning different values to the parameters.

Example 3 Optimization Considerations

A multi-resolution representation of an image I, referred to below as“the artist” image was produced according to some embodiments of thepresent invention using a Reduce operation featuring Gaussian weights,thereby forming a Gaussian pyramid. The process was applied as follows:B _(n)=Reduce(B _(n−1)),  (EQ. 16)where the finest resolution B₀ was set to be the original image:B₀=I.  (EQ. 17)

The obtained Gaussian pyramid for a set of six scaled images I₀, . . . ,I₅ is shown in FIGS. 13A-F.

A relative luminance q^(n) was then calculated according to therelation:q ^(n) =I ^(n) /I ^(n+1),  (EQ. 18)where I^(n) is the nth scaled image and I^(n+1) is an interpolatedversion of the n+1 scaled image obtained using the Expand operator:I^(n)=B_(n)I ^(n+1)=Expand(B _(n+1))  (EQ. 19)

The calculation of relative luminance resulted in a luminance pyramidcorresponding to the relative luminance levels q⁰, . . . , q⁴, shown inFIGS. 14A-E.

Once the luminance pyramid was obtained, a local contrast C^(k) wascalculated for each luminance level q. The local contrast was calculatedaccording to EQ. 6A above with ε=0.3. This resulted in a local contrastpyramid corresponding to the local contrast levels C⁰, . . . , C⁴, shownin FIGS. 15A-E.

Each of the local contrast levels C⁰, . . . , C⁴ was used forcalculating a image-specific exponent γ according to EQ. 5A, thusproviding a set of exponents γ⁰, . . . , γ⁴. The correspondingrepresentation pyramid is shown in FIGS. 16A-E. Thereafter, EQ. 4 wasemployed for calculating an effective saturation function R^(n) for eachluminance level q^(n), and each exponent γ^(n). This resulted in asaturation pyramid corresponding to the functions R⁰, . . . , R⁴, shownin FIGS. 17A-E.

The effective luminance levels were then combined according to EQ. 14B,to provide a combined image I_(combined). The original image I and thecombined image I_(combined) are shown in FIGS. 18A and 18B respectively.

Multi resolutions algorithms tend to suffer from halo artifacts thatappear around sharp edges. The present inventors contemplate bothembodiments in which a bilateral filter is used and alternativeembodiments in which a bilateral filter is not used. In the latterembodiments, considerations regarding the coarse contrast of the imageare preferably made. In various exemplary embodiments of the inventionthe processing technique applies one or more smoothing filter, such as,but not limited to, Gaussian filter and performs the enhancementadaptively, as further detailed hereinabove. For areas in the imagewhere the contrast is high, the contrast is not enhanced and may even bedecreases.

It was found by the present inventors that for images that arerelatively dark, the use of global gain may improve the quality of thefinal image. Thus, in some embodiments of the present invention, aglobal gain operation is applied in addition to the local contrastenhancement. This can be achieved in any technique known in the art. Insome embodiments, all of the image pixels are raised to the sameconstant power, which is optionally and preferably less than 1:I→I^(p)  (EQ. 20)

In various exemplary embodiments of the invention the gain exponent p isselected using an optimization procedure. For example, a set {p_(j)} ofcandidate values of the gain exponent can be selected and a score can beassigned to each candidate gain exponent p_(j) of the set. The gainexponent is then selected from the candidates responsively to theassigned scores. Typically, the gain exponent is the candidate havingthe highest score. Formally, denoting the score assigned to the jthcandidate gain exponent by S(I^(p) ^(j) ), the gain exponent can becalculated using the operation:

$\begin{matrix}{p = {\underset{p_{j}}{\arg\;\max}\left( {S\left( I^{p_{j}} \right)} \right)}} & \left( {{EQ}.\mspace{14mu} 21} \right)\end{matrix}$

A score suitable for the present embodiments is characteristic contrastof at least one of, more preferably at least some of, most preferablyall, the scaled images in the set. Such characteristic contrast isreferred to herein as a “Contrast Measure”.

In some embodiments, the Contrast Measure is a local-global contrast. Arepresentative example of a technique for calculating a local-globalcontrast suitable for the present embodiments is found in Rizzi et al.,“A Modified Algorithm for Perceived Contrast Measure in Digital Images,”CGIV 2008 and MCS'08 Final Program and Proceedings, the contents ofwhich are hereby incorporated by reference. For example, the ContrastMeasure can be calculated by summing the averages of the pyramid levels.

$\begin{matrix}{{{Contrast}\mspace{14mu}{Measure}\;(I)} = {\sum\limits_{n}\;{{Average}\left( {I_{n} - I_{n - 1}} \right)}}} & \left( {{EQ}.\mspace{14mu} 22} \right)\end{matrix}$

The summation in EQ. 22 is optionally and preferably performed over allthe elements in the set. EQ. 22 ensures that an image with high contrastyields higher values for the contrast measure, wherein an image with lowcontrast yields lower values for the contrast measure.

When the score S is enacted by the function Contrast Measure, EQ. 22 isapplied for each candidate p_(j), and the gain exponent p is preferablyselected using the operation:

$\begin{matrix}{p = {\underset{p_{j}}{\arg\;\max}{\left( {{Contrast}\mspace{14mu}{{Measure}\left( I^{p_{j}} \right)}} \right).}}} & \left( {{EQ}.\mspace{14mu} 23} \right)\end{matrix}$

The size of the set (namely, the number of scaled images in the setwhich is the number of different resolutions employed) can be selectedby the size of the image or by the amount of information in the image.

When the size of the set is selected based on the size of the image, itis optionally and preferably selected such that the coarsest resolutionis of a predetermined size, for example a 64×64 resolution or a 32×32resolution or a 16×16 resolution. Selecting the size of the set base onthe size of the image is particularly useful when a relatively smallportion of the input image is occupied by background.

When the size of the set is selected based on the amount of informationin the image, the size of the set is preferably determined during thebuildup of the set. In these embodiments, the amount of information thatis being added by each resolution is determined. Once the amount ofadded information is below a predetermined threshold, the set buildup isterminated. The amount of added information can be calculated, forexample, by counting the number of identifiable features (e.g., edges,distinguishable regions) in each scaled image. Selecting the size of theset base on the size of the image is particularly useful when arelatively large portion of the input image is occupied by background,and is generally preferred from the standpoint of computer resources.

Although the invention has been described in conjunction with specificembodiments thereof, it is evident that many alternatives, modificationsand variations will be apparent to those skilled in the art.Accordingly, it is intended to embrace all such alternatives,modifications and variations that fall within the spirit and broad scopeof the appended claims.

All publications, patents and patent applications mentioned in thisspecification are herein incorporated in their entirety by referenceinto the specification, to the same extent as if each individualpublication, patent or patent application was specifically andindividually indicated to be incorporated herein by reference. Inaddition, citation or identification of any reference in thisapplication shall not be construed as an admission that such referenceis available as prior art to the present invention. To the extent thatsection headings are used, they should not be construed as necessarilylimiting.

What is claimed is:
 1. A method of processing an image, comprising:obtaining an image decomposed into a set of scaled images, each beingcharacterized by a different image-scale; for each of at least somescaled image of said set, calculating a relative luminance between saidscaled image and another scaled image of said set, using intensities insaid scaled image and intensities in said another scaled image;processing each scaled image of said set using an adaptation procedurefeaturing an image-specific effective saturation function of saidrelative luminance, thereby providing a processed scaled image;combining at least some of the processed scaled images to provide acombined image; and outputting said combined image to a non-transitorycomputer readable medium.
 2. The method of claim 1, wherein saidobtaining comprising receiving the image and decomposing the image intosaid set of scaled images.
 3. The method of claim 1, wherein saiddecomposing comprises selecting a size of said set based on a size ofthe image.
 4. The method of claim 1, wherein said decomposing comprisesdetermining an amount of information in each scaled image being formed,and ceasing said decomposing when said amount of information is below apredetermined threshold.
 5. The method according to claim 1, wherein acharacteristic dynamic range of said combined image is lower than acharacteristic dynamic range of the original image.
 6. The methodaccording to claim 1, wherein said combining comprises multiplying. 7.The method according to claim 1, wherein said processing comprisesmodulating each relative luminance to provide a plurality of modulatedrelative luminance levels, wherein said combining comprises combiningsaid modulated relative luminance levels.
 8. The method according toclaim 1, wherein said set is an ordered set and wherein said relativeluminance is expressed as function of a ratio between said intensitiesin said scaled image and said intensities in said other scaled image. 9.The method according to claim 1, wherein said image-specific effectivesaturation function comprises an image-specific exponent, which is afunction of a local contrast within said scale-image.
 10. The methodaccording to claim 9, wherein said processing comprises modulating eachrelative luminance to provide a plurality of modulated relativeluminance levels, wherein said combining comprises combining saidmodulated relative luminance levels.
 11. The method according to claim10, wherein said modulating comprises selecting a relative luminancelevel such that two effective saturation functions corresponding todifferent image-specific exponent but the same scale be substantiallymatched.
 12. The method according to claim 9, wherein said localcontrast is calculated using a contrast-based adaptation procedureemployed for each picture-element of said scaled image.
 13. The methodaccording to claim 12, wherein said contrast-based adaptation procedurecalculates said local contrast based on a difference between a secondorder opponent receptive field function calculated for saidpicture-element and a second order opponent receptive field functioncalculated for nearby picture-elements.
 14. The method according toclaim 9, wherein said image-specific exponent is a decreasing functionof said local contrast.
 15. The method according to claim 14, whereinsaid image-specific exponent is a linear decreasing function of saidlocal contrast.
 16. The method according to claim 1, wherein saidimage-specific effective saturation function comprises a modulationfunction which is calculated based on a local contrast.
 17. The methodof claim 16, wherein said modulation function has higher values whensaid local contrast is low, and lower values when said local contrast ishigh.
 18. The method according to claim 1, further comprising employinga global gain operation for all scaled images of said set.
 19. Themethod according to claim 18, wherein said global gain operationfeatures a global gain exponent, and the method further comprisescalculating said global gain exponent using an optimization procedure.20. The method of claim 18, wherein said optimization procedurecomprises selecting a set of candidate gain exponents, assigning a scoreto each candidate gain exponent, and selecting said gain exponentresponsively to said score.
 21. The method of claim 20, wherein saidscore comprises a characteristic contrast.
 22. The method according toclaim 1, wherein said set is an ordered set and wherein said scaledimage and said other scaled image are adjacent images in said set. 23.The method according to claim 1, wherein the image is an HDR image. 24.The method according to claim 1, wherein the image is of at least onetype selected from the group consisting of a visible light image, astills image, a video image, an X-ray image, an infrared image, athermal image, a ultraviolet image, a computerized tomography (CT)image, a mammography image, a Roentgen image, a positron emissiontomography (PET) image, a magnetic resonance image, an ultrasoundimages, an impedance image, an elastography image, and a single photonemission computed tomography (SPECT) image.
 25. A method of capturingand displaying an image, comprising capturing an image of a scene andprocessing the image using the method according to claim
 1. 26. Themethod of claim 25, wherein said capturing the image comprises recordingradiation selected from the group consisting of visible light, infraredlight, ultraviolet light, X-ray radiation, radiofrequency radiation,microwave radiation and ultrasound radiation.
 27. A non-transitorycomputer software product, comprising a non-transitory computer-readablemedium in which program instructions are stored, which instructions,when read by a computer, cause the computer to execute the method ofclaim
 1. 28. A system for processing an image, the system comprises adata processor configured for: obtaining an image decomposed into a setof scaled images, each being characterized by a different image-scale;calculating, for each scaled image of said set, a relative luminancebetween said scaled image and another scaled image of said set usingintensities in said scaled image and intensities in said another scaledimage; processing each scaled image of said set using an adaptationprocedure featuring an image-specific effective saturation function ofsaid relative luminance, thereby providing a processed scaled image; andcombining at least some of the processed scaled images to provide acombined image.
 29. The system according to claim 28, wherein said dataprocessor is configured for receiving the image and decomposing theimage into said set of scaled images.
 30. An imaging system, comprisingan image capturing system and the system according to claim
 28. 31. Theimaging system of claim 30, wherein said capturing system is selectedfrom the group consisting of a digital camera, a video camera, a CMOSdigital camera, an infrared camera, a thermography device, an X-raycamera, a scanner, a microwave imaging, a computerized tomographyscanner, a single photon emission computed tomography device, a positronemission tomography device, a magnetic resonance imaging scanner, amammography scanner, an ultrasonic scanner, an impedance imaging system,an endoscopic imaging device, an elastography device, a radio telescope,a digital telescope, a digital microscope and a system for translatingan analog image to a digital image.